clc;clear all;close all;
global alpha betta depr sigmaC sigmaL phi govexp omega taumax N T ALPHA popweight
global kstst lstst cstst rstst wstst kstart kistart tauK0
global Delta1 Delta2 lamda transfer PVtaxes deficit kg lam kmax
global restrict deductible

%run mainbench.m;

load('benchN2.mat')

%find long-run steady state given initial guesses for lam, Delta1 and Delta2

kistart = kj;
kstart = sum(popweight.*kistart)+kg;
kmax=3;

%ALPHA=1;

lamda=lam;

restrict = 2
Delta1 = .3;
Delta2 = -Delta1;

%opts1=[1e-6 1 1 1e-5];
%kststguess = broydn('findstst',kstart,opts1)%,0,0)
kstst = fzero('findstst',[kstart-0.2 kmax],options)%,0,0)

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end

T=150; %we assume that the economy reaches the steady state after T periods

tauK0 = taumax;

%initial guesses 
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
    X(T+i) = lstst; %l
    X(2*T+i) = max(1-.05*i,0); %gamma
end;
deductible=-1;
X(3*T+1:3*T+1+2*(N-1)) = [Delta1 deductible lamda];

T=150; %we assume that the economy reaches the steady state after T periods
tauK0 = taumax;

ALPHall=[1 0.8 0.6 0.5 0.45 0.4 0.35 0.3 0.29 0.28 0.27];

save('solN2_ded.mat')

for ii=1:length(ALPHall)

ALPHA=ALPHall(ii);

if (ii>1)
    X=XX_ded(ii-1,:);
end

if ii==2 | ii==4 | ii==5 | ii==7 | ii==8
    X(3*T+3+N-2:3*T+3+2*(N-2))=0.6; %lambda
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
X1=X;

if ii==1
    X=X1;
    tolf0 = 1e-7;
    [X,check] = broydn('resid3v',X,tolf0)%,0,1);
    if size(X,1)>1
        X=X';
    end
    sum(resid3v(X).^2)

    options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
    [X,check] = fsolve('resid3v',X,options)
end

sumresid_all(ii)=sum(resid3v(X).^2)

XX_ded(ii,:)=X;
    
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = -Delta1
deductible = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
deductible_all(ii)=deductible;

kstst = fzero('findstst',[kstart-0.2 kmax],options)%,0,0)

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];
lnext = [X(T+2:2*T) lstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)- omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*(log(cistst) - omega*listst.^(1+sigmaL)/(1+sigmaL));
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC) - omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*((cistst.^(1-sigmaC)-1)/(1-sigmaC) - omega*listst.^(1+sigmaL)/(1+sigmaL));
end

if sum(V>=Vsq)==2
    PI(ii)=1
else
    PI(ii)=0
end
Vall_ded(ii,:)=V
Vsq

%tax rates
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
plot(tauKt)
[swi dur]=min((tauKt-taumax/2).^2);
duration(ii)=dur
deductible_all

save('solN2_ded.mat')

end

